Resistance and Robustness Visualizations and Statistics

Note on this code:

R version 3.5.1 (2018-07-02) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Sierra 10.12.6

Code written by Sara D. Williams

Resistance Results Visualization

The following code chunk is used to make Figure 3.

bleachmodel<-read.csv("Resistance.csv", header=T) #compiled results for Resistance from Bleaching model, see Metadata 2 for more information.
summary(bleachmodel)

       X                   Spatial.Scale       R               R_std        
 Min.   : 1.00   Caribbean        : 5    Min.   :0.09428   Min.   :0.00000  
 1st Qu.:18.25   Central Caribbean: 5    1st Qu.:0.19074   1st Qu.:0.01897  
 Median :35.50   Central Pacific  : 5    Median :0.25842   Median :0.02477  
 Mean   :35.50   Eastern Caribbean: 5    Mean   :0.34151   Mean   :0.02941  
 3rd Qu.:52.75   Eastern Pacific  : 5    3rd Qu.:0.50396   3rd Qu.:0.03397  
 Max.   :70.00   Global           : 5    Max.   :0.88822   Max.   :0.10159  
                 (Other)          :40                                       
                                 Simulation   Group       statlab  
 Network                              :14   Carib:15   A      : 7  
 Random Bipartite Degree Conserved    :14   Ind  :15   D      : 7  
 Random Bipartite Not-Degree Conserved:14   Main :20   C      : 6  
 Random Tolerances                    :14   Pac  :20   E      : 6  
 Shuffled Tolerances                  :14              I      : 6  
                                                       F      : 5  
                                                       (Other):33  

theme_set(theme_grey(base_size = 28)) 
#reorder simulations
bleachmodel$Simulation<-factor(bleachmodel$Simulation,levels=c("Network","Shuffled Tolerances","Random Tolerances","Random Bipartite Degree Conserved","Random Bipartite Not-Degree Conserved"))
#plot the global and ocean-basins
Main_oceans<-bleachmodel %>%
  filter(Group=='Main')
#reorder x axis
Main_oceans$Spatial.Scale<-factor(Main_oceans$Spatial.Scale,levels=c("Global","Pacific","Indian","Caribbean"))
main<-ggplot( Main_oceans, aes(x=as.factor(Spatial.Scale), y=R,fill=Simulation)) +
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=R-R_std, ymax=R+R_std), width=.2,position=position_dodge(.9))+
  scale_fill_brewer(palette="BuGn")+
  coord_cartesian(ylim=c(0,1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  xlab("Global and Oceans")+
  ylab("Resistance")+
  #labs(x="Global and Oceans",y="Resistance",size=10)+
  geom_text(aes( label=statlab,y=R+R_std+0.1),position = position_dodge(1),vjust =0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
#caribbean subregions
carib<-bleachmodel %>%
  filter(Group=='Carib')
c<-ggplot( carib, aes(x=as.factor(Spatial.Scale), y=R, fill=Simulation)) +
  geom_bar(position=position_dodge(), stat="identity",colour='black') +
  geom_errorbar(aes(ymin=R-R_std, ymax=R+R_std), width=.2,position=position_dodge(.9))+
  scale_fill_brewer(palette="BuGn")+coord_cartesian(ylim=c(0,1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="Caribbean Regions",y=" ")+
  geom_text(aes( label=statlab,y=R+R_std+0.1),position = position_dodge(1),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
#Pacific subregions
pac<-bleachmodel %>%
  filter(Group=='Pac')
p<-ggplot( pac, aes(x=as.factor(Spatial.Scale), y=R, fill=Simulation)) +
  geom_bar(position=position_dodge(), stat="identity", colour='black') +
  geom_errorbar(aes(ymin=R-R_std, ymax=R+R_std), width=.2,position=position_dodge(.9))+
  scale_fill_brewer(palette="BuGn")+coord_cartesian(ylim=c(0,1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="Pacific Regions",y=" ")+
  geom_text(aes( label=statlab,y=R+R_std+0.1),position = position_dodge(1),vjust=0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
#Indian subregions
ind<-bleachmodel %>%
  filter(Group=='Ind')
i<-ggplot( ind, aes(x=as.factor(Spatial.Scale), y=R, fill=Simulation)) +
  geom_bar(position=position_dodge(), stat="identity", colour='black') +
  geom_errorbar(aes(ymin=R-R_std, ymax=R+R_std), width=.2,position=position_dodge(.9))+
  scale_fill_brewer(palette="BuGn")+coord_cartesian(ylim=c(0,1))+theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="Indian Regions",y="Resistance")+
  geom_text(aes( label=statlab,y=R+R_std+0.1),position = position_dodge(1),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))

#Put them all in one plot
ggarrange(main, p, i, c, ncol=2, nrow=2, common.legend = TRUE, legend="bottom")

ggsave("figure3.jpg", plot = last_plot(), device = NULL, path = NULL,scale = 1, width = NA, height = NA, units = c("in"),dpi = 600)

Saving 24 x 16 in image

Robustness

The following code chunk is used to make Figure 4 E-H.

robustness<-read.csv("Robustness.csv") #results of robustness analyses, see metadata 2 for more info
summary(robustness)

    network         mean              std                   model          removed   
 C      : 30   Min.   :0.02041   Min.   :0.000000   bleach     : 42   both     : 56  
 cc     : 30   1st Qu.:0.53709   1st Qu.:0.008193   Degree_high: 42   hosts    : 56  
 cp     : 30   Median :0.64939   Median :0.027461   Degree_low : 42   links    :252  
 ec     : 30   Mean   :0.61124   Mean   :0.037324   LT_BH      : 42   symbionts: 56  
 ep     : 30   3rd Qu.:0.71870   3rd Qu.:0.054356   LT_BL      : 42                  
 G      : 30   Max.   :1.00000   Max.   :0.220333   LT_HL      : 42                  
 (Other):240                                        (Other)    :168                  
   type      R50_who       group2      group        statlab     connectance     
 link:252   both :140   carib :120   carib: 90          :366   Min.   :0.01000  
 node:168   hosts:140   Global: 30   ind  : 90   A      :  7   1st Qu.:0.03700  
            symbs:140   ind   :120   Main :120   B      :  7   Median :0.06650  
                        pac   :150   pac  :120   C      :  7   Mean   :0.07071  
                                                 D      :  6   3rd Qu.:0.08500  
                                                 E      :  6   Max.   :0.23600  
                                                 (Other): 21                    
     hosts         symbionts          links       
 Min.   : 14.0   Min.   : 13.00   Min.   :  43.0  
 1st Qu.: 31.0   1st Qu.: 29.00   1st Qu.:  84.0  
 Median : 83.5   Median : 40.50   Median : 209.5  
 Mean   :144.4   Mean   : 63.21   Mean   : 358.4  
 3rd Qu.:157.0   3rd Qu.: 74.00   3rd Qu.: 404.0  
 Max.   :685.0   Max.   :250.00   Max.   :1697.0  
                                                  

#get rid of the removal models that were tested but not used
data2<-robustness %>%
    filter(model!="LT_BH", model!="LT_HH",model!="LT_SH",model!="Tolerance_high")
  
Global_tot<-data2%>%
  filter(network=="G",removed!="hosts",removed!="symbionts",R50_who=="both")
Global_tot$removed<-factor(Global_tot$removed,levels=c("links"="links","nodes"="both"))
  
global<-ggplot( Global_tot, aes(x=as.factor(model), y=mean,fill=removed))
g_total<-global+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="Robustness, R50")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low","Degree_low","Degree_high"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Pacific_tot<-data2%>%
  filter(network=="P",removed!="hosts",removed!="symbionts",R50_who=="both")
Pacific_tot$removed<-factor(Pacific_tot$removed,levels=c("links"="links","nodes"="both"))
  
Pacific<-ggplot( Pacific_tot, aes(x=as.factor(model), y=mean,fill=removed))
p_total<-Pacific+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low","Degree_low","Degree_high"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Carib_tot<-data2%>%
  filter(network=="C",removed!="hosts",removed!="symbionts",R50_who=="both")
Carib_tot$removed<-factor(Carib_tot$removed,levels=c("links"="links","nodes"="both"))
Carib<-ggplot( Carib_tot, aes(x=as.factor(model), y=mean,fill=removed))
c_total<-Carib+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low","Degree_low","Degree_high"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Ind_tot<-data2%>%
  filter(network=="I",removed!="hosts",removed!="symbionts",R50_who=="both")
Ind_tot$removed<-factor(Carib_tot$removed,levels=c("links"="links","nodes"="both"))
Ind<-ggplot( Ind_tot, aes(x=as.factor(model), y=mean,fill=removed))
i_total<-Ind+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low","Degree_low","Degree_high"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
#Put them all in one plot
ggarrange(g_total,p_total,i_total,c_total, ncol=4, nrow=1, common.legend = FALSE, legend=NULL)

ggsave("figure4etof.jpg", plot = last_plot(), device = NULL, path = NULL,scale = 1, width = NA, height = NA, units = c("in"),dpi = 600)

Saving 24 x 16 in image

The following code chunk is used to make Appendix S2, Figure S1 D-F.

rr

myplot<-function(data2,N){ #data2=robustness results, N=network to plot
  data2<-data2 %>%
    filter(model!=\LT_BH\, model!=\LT_HH\,model!=\LT_SH\,model!=\Tolerance_high\)
  
  Global_tot<-data2%>%
    filter(network==N,removed!=\hosts\,removed!=\symbionts\,R50_who==\both\)
  Global_tot$removed<-factor(Global_tot$removed,levels=c(\links\,\both\,\hosts\,\symbionts\))
  
  global<-ggplot( Global_tot, aes(x=as.factor(model), y=mean,fill=removed))
  g_total<-global+
    geom_bar(position=position_dodge(),stat=\identity\,colour='black',mapping=aes(col=\red\)) +
    geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
    scale_fill_manual(values=c(\links\ = \#7fc97f\, \both\ = \#99d8c9\, \hosts\ = \#a6cee3\,\symbionts\=\#ffff99\), 
                      drop = FALSE)+
    
    #ylim(0,1)+
    scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
    theme(panel.background = element_blank())+ 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), axis.line = element_line(colour = \black\))+
    labs(x=\\,y=\Total Robustness\)+
    theme(axis.text.x = element_text(angle = 70, hjust = 1))+
    scale_x_discrete(limits=c(\Random_link\,\bleach\,\LT_BL\,\LT_HL\,\LT_SL\,\Random_node\,\Tolerance_low\,\Degree_low\,\Degree_high\),
                     labels=c(\Random_link\=\Random\,\bleach\=\Bleaching\,\Tolerance_low\=\Susceptible\,\Tolerance_high\=\High Tolerance\,\Degree_high\= \High Degree\,\Degree_low\=\Low Degree\,\Random_node\=\Random Node\,\LT_BH\=\High Avg. Link Tolerance\,\LT_BL\=\Susceptible\,\LT_HH\=\High Host Link Tolerance\,\LT_HL\=\Host Susceptible\,\LT_SH\=\High Symbiont Link Tolerance\,\LT_SL\=\Symbiont Susceptible\))+
    geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0)
  
  Global_h<-data2%>%
    filter(network==N,removed!=\hosts\,removed!=\both\,R50_who==\hosts\)
  Global_h$removed<-factor(Global_h$removed,levels=c(\links\,\both\,\hosts\,\symbionts\))
  global_h<-ggplot( Global_h, aes(x=as.factor(model), y=mean,fill=removed))
  g_hosts<-global_h+
    geom_bar(position=position_dodge(),stat=\identity\,colour='black',mapping=aes(col=\red\)) +
    geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
    scale_fill_manual(values=c(\links\ = \#7fc97f\, \both\ = \#99d8c9\, \hosts\ = \#a6cee3\,\symbionts\=\#ffff99\), 
                      drop = FALSE)+
    scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
    theme(panel.background = element_blank())+ 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), axis.line = element_line(colour = \black\))+
    labs(x=\\,y=\Host Robustness\)+
    theme(axis.text.x = element_text(angle = 70, hjust = 1))+
    scale_x_discrete(limits=c(\Random_link\,\bleach\,\LT_BL\,\LT_HL\,\LT_SL\,\Random_node\,\Tolerance_low\,\Degree_low\,\Degree_high\),
                     labels=c(\Random_link\=\Random\,\bleach\=\Bleaching\,\Tolerance_low\=\Susceptible\,\Tolerance_high\=\High Tolerance\,\Degree_high\= \High Degree\,\Degree_low\=\Low Degree\,\Random_node\=\Random Node\,\LT_BH\=\High Avg. Link Tolerance\,\LT_BL\=\Susceptible\,\LT_HH\=\High Host Link Tolerance\,\LT_HL\=\Host Susceptible\,\LT_SH\=\High Symbiont Link Tolerance\,\LT_SL\=\Symbiont Susceptible\))+
    geom_text(aes( label=statlab,y=mean+std+0.03),position = position_dodge(0.9),vjust =0)
  
  Global_s<-data2%>%
    filter(network==N,removed!=\symbionts\,removed!=\both\,R50_who==\symbs\)
  Global_s$removed<-factor(Global_s$removed,levels=c(\links\,\both\,\hosts\,\symbionts\))
  global_s<-ggplot( Global_s, aes(x=as.factor(model), y=mean,fill=removed))
  g_symbs<-global_s+
    geom_bar(position=position_dodge(),stat=\identity\,colour='black',mapping=aes(col=\red\)) +
    geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
    scale_fill_manual(values=c(\links\ = \#7fc97f\, \both\ = \#99d8c9\, \hosts\ = \#a6cee3\,\symbionts\=\#ffff99\), 
                      drop = FALSE)+
    scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
    theme(panel.background = element_blank())+ 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), axis.line = element_line(colour = \black\))+
    labs(x=\\,y=\Symbiont Robustness\)+
    theme(axis.text.x = element_text(angle = 70, hjust = 1))+
    scale_x_discrete(limits=c(\Random_link\,\bleach\,\LT_BL\,\LT_HL\,\LT_SL\,\Random_node\,\Tolerance_low\,\Degree_low\,\Degree_high\),
                     labels=c(\Random_link\=\Random\,\bleach\=\Bleaching\,\Tolerance_low\=\Susceptible\,\Tolerance_high\=\High Tolerance\,\Degree_high\= \High Degree\,\Degree_low\=\Low Degree\,\Random_node\=\Random Node\,\LT_BH\=\High Avg. Link Tolerance\,\LT_BL\=\Susceptible\,\LT_HH\=\High Host Link Tolerance\,\LT_HL\=\Host Susceptible\,\LT_SH\=\High Symbiont Link Tolerance\,\LT_SL\=\Symbiont Susceptible\))+
    geom_text(aes( label=statlab,y=mean+std+0.03),position = position_dodge(0.9),vjust =0)
  
  theme_set(theme_grey(base_size = 14)) 
  p<-ggarrange(g_total,g_hosts,g_symbs, ncol=3, nrow=1, common.legend = TRUE, legend=\bottom\)
  return(p)
}
myplot(robustness,\G\)

Randomization Tests

See Appendix S2 for description. Randomization tests are often called permutation tests, thus I use “permutation” within the code.

rr

mypermutation_twocomp<-function(data,netA,netB,levelA,levelB,nsims){
  #get data
  #levels(data$net_type)<-c(\hsrand_dtemp\, \net_dtemp\, \rand_dtemp\,\rbdc_dtemp\, \rbndc_dtemp\)
  #split data for two levels
  net1 <- data%>%
    filter(net_type==netA)
  net2 <- data%>%
    filter(net_type==netB)
  twocomps<-rbind(net1,net2)
  # initialize test
  combined_resistance <- c(net1$resistance,net2$resistance) #combines resistance values into a vector
  combined_labels <- c(net1$net_type,net2$net_type)#combines network type into a vector
  diff_obs <- mean(net1$resistance) - mean(net2$resistance)
  model <- anova(lm(twocomps$resistance ~ twocomps$net_type)) 
  obt.F <- model$\F value\[1]     # Our obtained F  statistic
  obt.p <- model$\Pr(>F)\[1]
  
  # permutation test
  diffs <- rep(NA, nsims)
  fstats<-rep(NA,nsims)
  for (i in 1:nsims) {
    #shuffle the labels
    shuffled_labels <- sample(combined_labels, replace = FALSE)
    twocomps$shuffled_labels<-shuffled_labels
    #shuffle the resistances
    shuffled_resistance<-sample(combined_resistance, replace = FALSE)
    twocomps$shuffled_resistance<-shuffled_resistance
    #get the differences in the means of the shuffled groups
    diffs[i] <- mean(shuffled_resistance[shuffled_labels == levelA]) -  mean(shuffled_resistance[shuffled_labels == levelB])
    #run a new anova
    newmodel <- anova(lm(twocomps$shuffled_resistance ~ twocomps$shuffled_labels)) 
    fstats[i] <- newmodel$\F value\[1]     # get a new F  statistic
  }
  #calculate the two-sided p-value
  # p-value = (number of more extreme differences than diff_obs)/nsims
  pval_diffs<-length(diffs[abs(diffs) >= abs(diff_obs)])/nsims
  pval_fstatdifs<-(length(fstats[abs(fstats) > abs(obt.F)])+1)/(1+nsims)
  return(c(diff_obs,obt.p,obt.F,pval_diffs,pval_fstatdifs))
}
resistance_perms<-function(file,nsims){
  #get data
  resist_data<-read.csv(file)
  #select columns
  selectdata<-resist_data%>%
    select(hsrand_dtemp,net_dtemp,rand_dtemp,rbdc_dtemp,rbndc_dtemp)
  #get into long format
  longdata <- gather(selectdata, key=net_type,value=resistance,factor_key = TRUE)
  #setup storage for permutation results
  resmat<-matrix(NA,nrow=0,ncol=7)
  pvals<-rep()
  #run permutations
  for (j in 1:(length(levels(longdata$net_type))-1)){
    netA<-levels(longdata$net_type)[j]
    for (i in (j+1):length(levels(longdata$net_type))){
      netB<-levels(longdata$net_type)[i]
      levelA<-j
      levelB<-i
      #print(c(netA,netB,levelA,levelB))
      permresults<-mypermutation_twocomp(longdata,netA,netB,levelA,levelB,1000)
      #print(results)
      #pvals<-append(pvals,permresults[5])
      permres<-c(permresults,netA,netB)
      resmat<-rbind(resmat,permres)
  }
}
  all_results<-cbind(resmat,p.adjust(as.numeric(resmat[,5]),\holm\)) #adjust p values using the holm correction
  results<-as.data.frame(all_results)
  colnames(results)<-c(\diff_obs\,\obt.p\,\obt.F\,\pval_diffs\,\pval_fstatdifs\,\net1\,\net2\,\p_adjust_holm\)
  return(results)
}
#TEST
resistance_perms(\TEST_resistance_perms.csv\,1000)

                     diff_obs                obt.p
permres   0.00412371134020573    0.702317494689414
permres.1  -0.436082474226805 6.33757407839271e-88
permres.2  -0.152577319587629 2.49559940772465e-38
permres.3  -0.121649484536083 5.37670669053595e-26
permres.4  -0.440206185567011 3.32885050387969e-86
permres.5  -0.156701030927835 5.16592722639757e-37
permres.6  -0.125773195876289 1.85944491914943e-25
permres.7   0.283505154639176 2.96101092082189e-63
permres.8   0.314432989690722 3.42201761802395e-67
permres.9  0.0309278350515463 0.000594160717217541
                      obt.F pval_diffs
permres   0.146508966043483      0.776
permres.1   1318.0773480663          0
permres.2  268.554789272025          0
permres.3  151.108297535609          0
permres.4  1257.08014938236          0
permres.5  254.297638156382          0
permres.6  146.730745532962          0
permres.7  644.188110026627          0
permres.8  726.876119160036          0
permres.9  12.1964991530209      0.001
                pval_fstatdifs         net1
permres      0.725274725274725 hsrand_dtemp
permres.1 0.000999000999000999 hsrand_dtemp
permres.2 0.000999000999000999 hsrand_dtemp
permres.3 0.000999000999000999 hsrand_dtemp
permres.4 0.000999000999000999    net_dtemp
permres.5 0.000999000999000999    net_dtemp
permres.6 0.000999000999000999    net_dtemp
permres.7 0.000999000999000999   rand_dtemp
permres.8 0.000999000999000999   rand_dtemp
permres.9    0.001998001998002   rbdc_dtemp
                 net2       p_adjust_holm
permres     net_dtemp   0.725274725274725
permres.1  rand_dtemp 0.00999000999000999
permres.2  rbdc_dtemp 0.00999000999000999
permres.3 rbndc_dtemp 0.00999000999000999
permres.4  rand_dtemp 0.00999000999000999
permres.5  rbdc_dtemp 0.00999000999000999
permres.6 rbndc_dtemp 0.00999000999000999
permres.7  rbdc_dtemp 0.00999000999000999
permres.8 rbndc_dtemp 0.00999000999000999
permres.9 rbndc_dtemp 0.00999000999000999

Robustness vs. Connectance

See Appendix S2 for description.

rr

#load data
robustness<-read.csv(\Robustness.csv\)
#check distribution of connectance
plot(density((log(robustness$connectance)))) #It's wonky.

rr

#make it easy to plot removal model of choice vs connectance with a linear fit
my_RClinreg_plot<-function(choice){
  data<-robustness %>%
    filter(model!=\LT_BH\, model!=\LT_HH\,model!=\LT_SH\,model!=\Tolerance_high\,R50_who==\both\,model==choice)
  p<-ggplotRegression(lm(mean ~ connectance, data = data),choice)
  p+geom_point(y=data$mean,x=data$connectance,aes(color=data$group),size=6)
  
}
#run a Kendall's coefficient of rank correlation test on removal model of choice
my_cortest<-function(choice){
  data<-robustness %>%
    filter(model!=\LT_BH\, model!=\LT_HH\,model!=\LT_SH\,model!=\Tolerance_high\,R50_who==\both\,model==choice)
  ct<-cor.test(data$mean,data$connectance,method=\kendall\)
  return(ct)
}
#run for different removal models 
my_RClinreg_plot(\bleach\)

rr

my_cortest(\bleach\) 

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 67, p-value = 0.01928
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.4725275 

rr

my_RClinreg_plot(\Random_link\)

rr

my_cortest(\Random_link\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 52, p-value = 0.5183
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.1428571 

rr

my_RClinreg_plot(\LT_BL\)

rr

my_cortest(\LT_BL\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 55, p-value = 0.3308
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.2087912 

rr

my_RClinreg_plot(\LT_HL\)

rr

my_cortest(\LT_HL\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 56, p-value = 0.2792
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.2307692 

rr

my_RClinreg_plot(\LT_SL\)

rr

my_cortest(\LT_SL\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 42, p-value = 0.7472
alternative hypothesis: true tau is not equal to 0
sample estimates:
        tau 
-0.07692308 

rr

my_RClinreg_plot(\Degree_high\)

rr

my_cortest(\Degree_high\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 72, p-value = 0.003025
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.5824176 

rr

my_RClinreg_plot(\Degree_low\)

rr

my_cortest(\Degree_low\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 27, p-value = 0.04718
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.4065934 

rr

my_RClinreg_plot(\Random_node\)

rr

my_cortest(\Random_node\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 47, p-value = 0.9145
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
0.03296703 

rr

my_RClinreg_plot(\Tolerance_low\)

rr

my_cortest(\Tolerance_low\)

    Kendall's rank correlation tau

data:  data$mean and data$connectance
T = 56, p-value = 0.2792
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.2307692 

Symbiont node thermal tolerances: picking missing tolerances

Code for Appendix S1, Figure S1: Frequency distribution of Symbiodinium thermal tolerance scores adapted from Swain et al. (2017); distribution is colored by tolerance range, red is highly susceptible, orange is medium tolerance, and yellow is high tolerance. Swain et al. (2017) provides a framework for a consensus of Symbiodinium thermotolerance ranks developed from rank-aggregation methods. Their ranking scheme orders Symbiodinium phylotypes from 0-100, but the rank values are not indicative of total magnitude differences in thermotolerance. To determine tolerances of the unlisted symbiont types in our network, rank values were randomly drawn from the high, medium, and low thermal tolerance frequency distributions in the relative proportions of the clades represented in those distributions. Thus, for each simulation of either the bleaching or different removal models described below, the symbiont tolerances varied within a set distribution (Figure 1D).

rr

### Exploring data from Swain et al 2016a ####
data<-read.csv(\SwainSymbsR.csv\)
#frequency distribution --> inset from figure 1 from paper
bins=seq(0,100,by=2)
#split the frequency distribution into 3 subsets defined by swain et al 16
low <- data[which(data$ScoreR_k<=17),]
high <- data[which(data$ScoreR_k>=33),]
middle<- data[which(data$ScoreR_k>17 & data$ScoreR_k<33),]
#High hist
hist(high$ScoreR_k, breaks=bins, col=\yellow\, ylim=c(0,12), xlab=\thermotolerance score\, main=\Frequency Distributions of Thermotolerance Scores\)
#middle hist
hist(middle$ScoreR_k, breaks=bins, col=\orange\, add=T)

#low hist
hist(low$ScoreR_k, breaks=bins, col=\red\, add=T)

Code for fitting the above thermal toelrance distributions distributions

rr

### Fit the Thermotolerance Distributions ####
#first need the mixtools package
if(!(\mixtools\ %in% installed.packages())){install.packages(\mixtools\)}

Warning in install.packages :
  cannot open URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.5/PACKAGES.rds': HTTP status was '404 Not Found'
trying URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.5/mixtools_1.1.0.tgz'
Content type 'application/x-gzip' length 1413097 bytes (1.3 MB)
==================================================
downloaded 1.3 MB

The downloaded binary packages are in
    /var/folders/9_/h1rwfkm916nc6pcl1dsx6r6h0000gn/T//Rtmpfw2FHQ/downloaded_packages

rr

library(\mixtools\)

mixtools package, version 1.1.0, Released 2017-03-10
This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.

rr

#Use the expectation Maximization (EM) Approach
em <- normalmixEM(data$ScoreR_k, arbvar = TRUE,epsilon = 1e-03, k = 3)

number of iterations= 8 

rr

# Estimated means
em$mu

[1]  7.993676 34.266710 54.339110

rr

# Estimated standard deviations
em$sigma

[1]  2.531190 17.196483  1.094525

rr

###RESULTS
#> em$mu
#[1]  7.983029 29.408570 51.900339
#> em$sigma
#[1]  2.711480  9.500587 18.006626
#do 10000 simulations ####
runs <- 10000
sims.low <- rnorm(runs,mean=7.983029,sd=2.711480)
sims.middle <-rnorm(runs,mean=29.408570,sd=9.500587)
sims.high<-rnorm(runs,mean=51.900339,sd=18.006626)
#High hist
hist(sims.high, col=\yellow\)
#middle hist
hist(sims.middle, col=\orange\, add=T)

#low hist
hist(sims.low, col=\red\, add=T)

Code for setting up for toelrance function

rr

#VALUES THAT ARE NEEDED #### 
##These are just the numbers of symbs from each clade in each freq level from Swain
swain_A<-8
swain_B<-8
swain_C<-71
swain_D<-19
swain_E<-1
swain_F<-3
low_A<-1
middle_A<-3
high_A<-4
prob_low_A<-low_A/swain_A
prob_middle_A<-middle_A/swain_A
prob_high_A<-high_A/swain_A
low_B<-3
middle_B<-4
high_B<-1
prob_low_B<-low_B/swain_B
prob_middle_B<-middle_B/swain_B
prob_high_B<-high_B/swain_B
low_C<-24
middle_C<-19
high_C<-28
prob_low_C<-low_C/swain_C
prob_middle_C<-middle_C/swain_C
prob_high_C<-high_C/swain_C
low_D<-3
middle_D<-4
high_D<-12
prob_low_D<-low_D/swain_D
prob_middle_D<-middle_D/swain_D
prob_high_D<-high_D/swain_D
low_E<-1
prob_low_E<-1
prob_middle_E<-0
prob_high_E<-0
high_F<-3
prob_low_F<-0
prob_middle_F<-0
prob_high_F<-1
#OK so these are the number of symbs in each clade that dont have tolerances 
net_A<-6
net_B<-13
net_C<-160
net_D<-6
net_E<-0
net_F<-0

Function for getting the tolerances

rr

get_tols<-function(size,lowprob,midprob,highprob,runs,sims.low,sims.middle,sims.high){
  #make a sample distribution where 1 corresponds to the low thermaltolerance,
  # 2 is the medium and 3 is high. 1:3 are put into the sample distribution 
  #based on the probabilities in the function input that are clade specific
  sampledistr<-sample(x = 1:3, size = size, prob = c(lowprob, midprob, highprob), replace=TRUE)
  tolerance<-numeric(length=size)
  for (i in 1:size) { #for each number in the sampledistr do:
    if (sampledistr[i]==1){
      tolerance[i]<-sample(x=sims.low,size=1) #pull from low
    }
    else if (sampledistr[i]==2){ #pull from middle
      tolerance[i]<-sample(x=sims.middle,size=1)
    }
    else 
      tolerance[i]<-sample(x=sims.high,size=1) #pull from high
  }
  tolerance<-sqrt(tolerance)/10 #take the sqrt because swain's tolerances are not indicative of magnitude and then divide by 10 to get a number from 0-1
  return(tolerance)
}

Make 100 simulations worth of tolerance files

rr

trials=1 #actually did 100, but don't want to make a ton of new files, so here is 1 example file.
lst=seq(1:trials)
for(i in seq_along(lst)){
  Ctols<-get_tols(net_C,prob_low_C,prob_middle_C,prob_high_C,10000,sims.low,sims.middle,sims.high)
  Atols<-get_tols(net_A,prob_low_A,prob_middle_A,prob_high_A,10000,sims.low,sims.middle,sims.high)
  Btols<-get_tols(net_B,prob_low_B,prob_middle_B,prob_high_B,10000,sims.low,sims.middle,sims.high)
  Dtols<-get_tols(net_D,prob_low_D,prob_middle_D,prob_high_D,10000,sims.low,sims.middle,sims.high)
  Ctols<-as.data.frame(Ctols)
  Atols<-as.data.frame(Atols)
  Btols<-as.data.frame(Btols)
  Dtols<-as.data.frame(Dtols)
  colnames(Ctols) <- c(\tolerance\)
  colnames(Btols) <- c(\tolerance\)
  colnames(Dtols) <- c(\tolerance\)
  colnames(Atols) <- c(\tolerance\)
  newscores<-rbind(Atols,Btols,Ctols,Dtols)
  name<-paste('trial',i,\.csv\,sep=\\)
  
  write.csv(newscores,name,row.names=FALSE)
}

Now for the shuffled tolerances:

rr

# Now get the toleranes for the Shuffled Symbionts/ Shuffled Random null model ####
#for the random set of tolerances where they are all pulled equally from the 3 distributions
runs=10000
sims.low <- rnorm(runs,mean=7.983029,sd=2.711480)
sims.low<-subset(sims.low,sims.low>=0)
sims.middle <-rnorm(runs,mean=29.408570,sd=9.500587)
sims.middle<-subset(sims.middle,sims.middle>=0)
sims.high<-rnorm(runs,mean=51.900339,sd=18.006626) 
sims.high<-subset(sims.high,sims.high<=100)
sims.high<-subset(sims.high,sims.high>=0)
#so the following function does what the old one did but also specifies
#if using probabilies determined from data or if pulling all values from 
#the distributions with the same probability (1/3)
allthedata<-function(probability,trials,sims.low,sims.middle,sims.high,net_A,net_B,net_C,net_D){
  if (probability==1){
    swain_A<-8
    swain_B<-8
    swain_C<-71
    swain_D<-19
    swain_E<-1
    swain_F<-3
    
    low_A<-1
    middle_A<-3
    high_A<-4
    prob_low_A<-low_A/swain_A
    prob_middle_A<-middle_A/swain_A
    prob_high_A<-high_A/swain_A
    
    low_B<-3
    middle_B<-4
    high_B<-1
    prob_low_B<-low_B/swain_B
    prob_middle_B<-middle_B/swain_B
    prob_high_B<-high_B/swain_B
    
    low_C<-24
    middle_C<-19
    high_C<-28
    prob_low_C<-low_C/swain_C
    prob_middle_C<-middle_C/swain_C
    prob_high_C<-high_C/swain_C
    
    low_D<-3
    middle_D<-4
    high_D<-12
    prob_low_D<-low_D/swain_D
    prob_middle_D<-middle_D/swain_D
    prob_high_D<-high_D/swain_D
    
    low_E<-1
    prob_low_E<-1
    prob_middle_E<-0
    prob_high_E<-0
    
    high_F<-3
    prob_low_F<-0
    prob_middle_F<-0
    prob_high_F<-1
    
  }
  else
    prob_low_C=prob_middle_C=prob_high_C=prob_low_A=prob_middle_A=prob_high_A=prob_low_B=prob_middle_B=prob_high_B=prob_low_D=prob_middle_D=prob_high_D=(1/3)
  
  #net_A<-8
  #net_B<-18
  #net_C<-177
  #net_D<-7
  #net_E<-0
  #net_F<-0
  lst=seq(1:trials)
  for(i in seq_along(lst)){
    Ctols<-get_tols(net_C,prob_low_C,prob_middle_C,prob_high_C,10000,sims.low,sims.middle,sims.high)
    Atols<-get_tols(net_A,prob_low_A,prob_middle_A,prob_high_A,10000,sims.low,sims.middle,sims.high)
    Btols<-get_tols(net_B,prob_low_B,prob_middle_B,prob_high_B,10000,sims.low,sims.middle,sims.high)
    Dtols<-get_tols(net_D,prob_low_D,prob_middle_D,prob_high_D,10000,sims.low,sims.middle,sims.high)
    Ctols<-as.data.frame(Ctols)
    Atols<-as.data.frame(Atols)
    Btols<-as.data.frame(Btols)
    Dtols<-as.data.frame(Dtols)
    colnames(Ctols) <- c(\tolerance\)
    colnames(Btols) <- c(\tolerance\)
    colnames(Dtols) <- c(\tolerance\)
    colnames(Atols) <- c(\tolerance\)
    newscores<-rbind(Atols,Btols,Ctols,Dtols)
    name<-paste('rtrial',i,\.csv\,sep=\\)
    
    write.csv(newscores,name,row.names=FALSE)
  }
}
#size of the clades in the network
net_A<-10
net_B<-20
net_C<-211
net_D<-9
allthedata(0,1,sims.low,sims.middle,sims.high,net_A,net_B,net_C,net_D)

Revised Figures

Figure 4, now without the degree removal models

#,fig.height=2,fig.width=4.5}
robustness<-read.csv("Robustness.csv") #results of robustness analyses, see metadata 2 for more info
summary(robustness)

    network         mean              std                   model          removed   
 C      : 30   Min.   :0.02041   Min.   :0.000000   bleach     : 42   both     : 56  
 cc     : 30   1st Qu.:0.53709   1st Qu.:0.008193   Degree_high: 42   hosts    : 56  
 cp     : 30   Median :0.64939   Median :0.027461   Degree_low : 42   links    :252  
 ec     : 30   Mean   :0.61124   Mean   :0.037324   LT_BH      : 42   symbionts: 56  
 ep     : 30   3rd Qu.:0.71870   3rd Qu.:0.054356   LT_BL      : 42                  
 G      : 30   Max.   :1.00000   Max.   :0.220333   LT_HL      : 42                  
 (Other):240                                        (Other)    :168                  
   type      R50_who       group2      group        statlab     connectance     
 link:252   both :140   carib :120   carib: 90          :366   Min.   :0.01000  
 node:168   hosts:140   Global: 30   ind  : 90   A      :  7   1st Qu.:0.03700  
            symbs:140   ind   :120   Main :120   B      :  7   Median :0.06650  
                        pac   :150   pac  :120   C      :  7   Mean   :0.07071  
                                                 D      :  6   3rd Qu.:0.08500  
                                                 E      :  6   Max.   :0.23600  
                                                 (Other): 21                    
     hosts         symbionts          links       
 Min.   : 14.0   Min.   : 13.00   Min.   :  43.0  
 1st Qu.: 31.0   1st Qu.: 29.00   1st Qu.:  84.0  
 Median : 83.5   Median : 40.50   Median : 209.5  
 Mean   :144.4   Mean   : 63.21   Mean   : 358.4  
 3rd Qu.:157.0   3rd Qu.: 74.00   3rd Qu.: 404.0  
 Max.   :685.0   Max.   :250.00   Max.   :1697.0  
                                                  

#get rid of the removal models that were tested but not used
data2<-robustness %>%
    filter(model!="LT_BH", model!="LT_HH",model!="LT_SH",model!="Tolerance_high",model!="Degree_high",model!="Degree_low")%>%
  droplevels()
levels(data2$model)

[1] "bleach"        "LT_BL"         "LT_HL"         "LT_SL"         "Random_link"  
[6] "Random_node"   "Tolerance_low"

Global_tot<-data2%>%
  filter(network=="G",removed!="hosts",removed!="symbionts",R50_who=="both")
Global_tot$removed<-factor(Global_tot$removed,levels=c("links"="links","nodes"="both"))
  
global<-ggplot( Global_tot, aes(x=as.factor(model), y=mean,fill=removed))
g_total<-global+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="Robustness, R50")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Pacific_tot<-data2%>%
  filter(network=="P",removed!="hosts",removed!="symbionts",R50_who=="both")
Pacific_tot$removed<-factor(Pacific_tot$removed,levels=c("links"="links","nodes"="both"))
  
Pacific<-ggplot( Pacific_tot, aes(x=as.factor(model), y=mean,fill=removed))
p_total<-Pacific+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Carib_tot<-data2%>%
  filter(network=="C",removed!="hosts",removed!="symbionts",R50_who=="both")
Carib_tot$removed<-factor(Carib_tot$removed,levels=c("links"="links","nodes"="both"))
Carib<-ggplot( Carib_tot, aes(x=as.factor(model), y=mean,fill=removed))
c_total<-Carib+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
Ind_tot<-data2%>%
  filter(network=="I",removed!="hosts",removed!="symbionts",R50_who=="both")
Ind_tot$removed<-factor(Carib_tot$removed,levels=c("links"="links","nodes"="both"))
Ind<-ggplot( Ind_tot, aes(x=as.factor(model), y=mean,fill=removed))
i_total<-Ind+
  geom_bar(position=position_dodge(),stat="identity",colour='black',mapping=aes(col="red")) +
  geom_errorbar(aes(ymin=mean-std, ymax=mean+std), width=.2,position=position_dodge(.9))+
  scale_fill_manual(values=c("links" = "#7fc97f", "both" = "#99d8c9", "hosts" = "#a6cee3","symbionts"="#ffff99"), drop = FALSE)+
  #ylim(0,1)+
  scale_y_continuous(limits=c(0,1.05),expand = c(0,0),breaks=seq(0, 1, 0.1))+
  theme(panel.background = element_blank())+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(x="",y="")+
  theme(legend.position="none",axis.text.x = element_text(angle = 70, hjust = 1))+
  scale_x_discrete(limits=c("Random_link","bleach","LT_BL","LT_HL","LT_SL","Random_node","Tolerance_low"),labels=c("Random_link"="Random","bleach"="Bleaching","Tolerance_low"="Susceptible","Tolerance_high"="High Tolerance","Degree_high"= "High Degree","Degree_low"="Low Degree","Random_node"="Random Node","LT_BH"="High Avg. Link Tolerance","LT_BL"="Susceptible","LT_HH"="High Host Link Tolerance","LT_HL"="Host Susceptible","LT_SH"="High Symbiont Link Tolerance","LT_SL"="Symbiont Susceptible"))+
  geom_text(aes( label=statlab,y=mean+std+0.05),position = position_dodge(0.9),vjust = 0,size=8)+theme(
  axis.title.x = element_text(size = 24),
  axis.text.x = element_text(size = 20),
  axis.title.y = element_text(size = 24),
  axis.text.y=element_text(size=20))
#Put them all in one plot
ggarrange(g_total,p_total,i_total,c_total, ncol=4, nrow=1, common.legend = FALSE, legend=NULL)

ggsave("figure4etof.jpg", plot = last_plot(), device = NULL, path = NULL,scale = 1, width = NA, height = NA, units = c("in"),dpi = 600)

Saving 24 x 16 in image